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The aim of this paper is to present a kinetic numerical scheme for the computa- 
tions of transient pressurised flows in closed water pipe with non uniform sections. 
Firstly, we detail the derivation of the mathematical model in curvilinear coordi- 
nates and we performe a formal asymptotic analysis. The obtained system is written 
as a conservative hyperbolic partial differential system of equations. We obtain a 
■ kinetic interpretation of this system and we build the corresponding kinetic scheme 

, based on an upwinding of the source terms written as the gradient of a "pseudo 

altitude". The validation is lastly performed in the case of a water hammer in an 
uniform pipe: we compare the numerical results provided by an industrial code 
^ . used at EDF-CIH (France), which solves the Allievi equations (the commonly used 

<^ [ equation for pressurised flows in pipe) by the method of characteristics, with those 

^ ' of the kinetic scheme. To validate the contracting or expanding cases, we compare 

the presented technique to the equivalent pipe method in the case of an immediate 
^ i flow shut down in a quasi-frictionless cone-shaped pipe. 
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X : 1 Introduction 



The presented work takes place in a more general project: the modelization of 
unsteady mixed flows in any kind of closed domain taking into account the cavitation 
problem and air entrapment. We are interested in flows occuring in closed pipe of 
non uniform sections, where some parts of the flow can be free surface (it means 
that only a part of the pipe is filled) and other parts are pressurised (it means 
that the pipe is full-filled). The transition phenomenon, between the two types 
of flows, occurs in many situation such as storm sewers, waste or supply pipes in 
hydroelectric installation. It can be induced by sudden change in the boundary 
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conditions as failure pumping. During this process, the pressure can reach severe 
values and cause damages. 

The classical Shallow Water equations are commonly used to describe free sur- 
face flows in open channel. They are also used in the study of mixed flows using 
the Preissman slot artefact (see for example [IKII]). However, this technic does not 
take into account depressurisation phenomenon which occurs during a water ham- 
mer. We can also cite the Allievi equations which are commonly used to describe 
pressurised flows. Nonetheless, the non conservative form is not well adapted to a 
natural coupling with the Shallow Water equations (contrary to the one presented 
in m). 

The model for the unsteady mixed water flows in closed water pipes and a finite 
volume discretization has been previously studied by two of the authors (S] and a 
kinetic formulation has been proposed in [6]. This paper tends to extend naturally 
the work in [B] in the case of closed pipes with non uniform sections. 

We establish, in Section [2], the model for pressurised flows in curvilinear coordi- 
nates and recall some classical properties of this model. Rewritting the source terms 
due to both topography and geometry into a single one that we called pseudo-altitude 
term, we get a model close to the presented one by the authors in [9]. In Section [3l 
we present the kinetic formulation of this model that will be useful to show the main 
properties of the numerical scheme. The last part is devoted to the construction of 
the kinetic scheme: the upwinding of the source term due to the pseudo topography 
is performed in a close manner described by Perthame et al. |9] using an energetic 
balance at microscopic level. We have used the generalized characteristics method 
to extend the works in [6] to the kinetic scheme with pseudo-reflections. 

Finally, we present in Section a numerical validation of this study in the uni- 
form case by the comparison between the resolution of this model and the resolution 
of the Allievi equation solved by the industrial code belier used at Center in Hy- 
draulics Engineering of Electricite De France (EDF) [T2j for the case of critical water 
hammer tests. The validation in non uniform pipes is performed in the case of an 
immediate flow shut down in a quasi-frictionless cone-shaped pipe. The results are 
compared to the equivalent pipe method [I]. 



2 Formal Derivation of the model 

The presented model is derived from the 3D compressible Euler system written in 
curvilinear coordinates, then integrated over sections orthogonal to the main flow 
axis (see below). We neglect the second and third equation of the conservation of 
the momentum and we get an unidirectionnal model. Then, an asymptotic analysis 
is performed to get a model close to the Shallow Water model (to a future coupling 
for the study of unsteady mixed flows [1]). 

2.1 The Euler system in curvilinear coordinates 

The 3D Euler system in the cartesian coordinates is written as follows 

dtp + div(p[7) = 0, (1) 

dt{pU) + d\w{pU (^U) + Vp = F, (2) 
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where U {t, x,y, z) and p{t,x,y, z)) denotes the velocity with components {u,v,w) 
and the density respectively. p{t, x, y, z) is the scalar pressure and F the exterior 
strenght of gravity. 

We define the domain of the flow as the union of sections Vl{x) (assumed to be 
simply connected compact sets) orthogonal to some plane curve with parametriza- 
tion (x,0,6(x)) in a convenient cartesian reference frame (O, i , j , k) where k 
follows the vertical direction; b{x) is then the elevation of the point u;{x,0,b{x)) 
over the plane (O, i , j ) (see FiG. [T]). The curve may be, for instance, the axis 
spanned by the center of mass of each orthogonal section Q{x) to the main mean 
flow axis, especially in the case of a piecewise cone-shaped pipe. Notice that we 
consider only the case of rigid pipe: the sections are only x-dependent. 
To see the local effect induced by the geometry due to the change of sections and/or 
slope, we write the 3D compressible Euler system in the curvilinear coordinates. To 
this end, let us introduce the curvilinear variable defined by 



we denote by Z the altitude of any fluid particle M in the Serret-Prenet basis 
{T , N, B) at point ll'(x, 0, b{x)): T is the tangent vector, the normal vector and 
B the binormal vector (see FiG. Then we perform the following transformation 
T : (x, y, z) {X, Y, Z) and we use the following lemma (whose proof can be found 




where xq is an arbitrary abscissa. We set y = Y and 




T 



Then, for any vector field $ one has, 



In particular, for any scalar function f , one has 
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Figure 1: Geometric characteristics of the pipe 

Let ([/, V, WY be the components of the velocity vector in the {X, Y, Z) coordinates 
in such a way that the flow is orthogonal to the sections Q{x). Let R be the matrix 



defined by R 
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cos 



then: 




Applying Lemma [2.11 to the mass conservation equation, we get 

Jidtp + divipU)) = 



where 



dtiJp) + dxipU) + dripJV) + dzipJW) = 



J =det 



i-z-^e] cose sine ^ 

dX 



\ 





d 



1 







1- Z—e]sine cos6' 
dX 



(3) 



(4) 



To get the unidirectionnal Shallow Water-like equations, we suppose that the mean 
flow follows the X-axis. Hence, we neglect the second and third equation for the 
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conservation of the momentum. Therefore, we only perform the curvilinear trans- 
formation for the first conservation equation. To this end, multiplying the conser- 

cos ( 

vation of the momentum equation of System ([2]) by J I I and using Lemma 




O yields: 




dtipU) + div(pC/ (g)U)+Vp = -pVCg.OM) 

It may be rewritten as: 

dtiJpU) + dxipU^) + dripJUV^) + dzipJUW) + dxp 

= -pjg sine + pUW— (cose) 
dX 



where OM denotes the position of any particule M in the local basis 
(t,N,13) at point uj{x,0,b{x)). 

Finally, in the {X, Y, Z) coordinates the system reads: 

dt{jp) + dx{pU) + dYipJV) + dzipJW) = 

dtiJpU) + dxipU^) + dripJUV^) + dz{pJUW) + dxP (6) 

= -pJg sin e + pUW -^{cos e) 
dX 

Remark 2.1 Notice that k{X) = -p^e is the algebric curvature of the axis at uj{x) 

dX 

and the function J{X, Y, Z) = 1 — Zk{X) only depends on variables X, Z. Morever, 
we assume J > d m. Vip which corresponds to a reasonnable geometric hypothesis. 
Consequently, J defines a diffeomorphism and thus the performed transformation 
is admissible. 



We recall that the main objective is to obtain a formulation close to the Shallow 
Water equation in order to couple the two models in a natural way (in a close 
manner described in jSj). The direct integration of Equations ^ over Vt,{x) gives a 
model which is not useful, due to the term J, to perform a natural coupling with the 
Shallow Water model [1] for non uniform pipes. Setting e = H/L a small parameter 
(where H and L are two characteristics dimensions along k and i axis respec- 
tively), we get J = 1 + 0{e). We also assume that the characteristic dimension 
along the j axis is the same as k . We introduce the others characteristics dimen- 
sions T, P, [/, V , W for time, pressure and velocity repectively and the dimensionless 
quantities as follows: 

U = U/U, V = eV/U, W = eW/U, 

X = X/L, Y = Y/H, z = Z/H, p = p/p, e = e,p = p. 

In the sequel, we set P = and L = TU (i.e. we consider only laminar flow). 
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Under these hypothesis J{X,Y,Z) = J{X,Y,Z) = 1 — eZ—:^9. So, the rescaled 

dX 

system ([6]) reads: 



d^jp) + d^ipU) + dyiJpV) + d^iJpW) = 

d^iJUp) + d^ipU^) + dyiJpUV) + d^iJpUW) + %p 

,sin6' Zdj^{cos9) 



(7) 



epUWp{X)-p- 2 



with F, 



U 



/P7 



the Proude number along the i axis and the k or j axis 



where M is any generic variable. 

Formally, when e vanishes the system reduces to: 



dr{p} + d^{pU)+dy{pV) + d^{pW) = 



_sm y 
Z5^(cos0) 



(8) 



Finally, the system in variables (X, y, Z) that describes the slope variation and the 
section variation in a closed pipe reads: 



dt{p) + dx{pU) + dY{pV) + dz{pW) = 

dt{Up) + dx{pU^) + dripUV) + dz{pUW) + dxp = 



-pg sm 



(9) 



Remark 2.2 To take into account the friction, we add the source term —pgSfT 
(described above) in the momentum equation. 



2.2 Shallow Water-like equations in closed pipe 

In the following, we use the linearized pressure law p = Pa + ^ ^ (see e.g. [TTlfTS] ) 

PPo 

in which po represents the density of the fluid at atmospheric pressure pa and /? the 
water compressibility coefficient equal to 5.0 10"^'' m^.A^~^ in practice. The sonic 
speed is then given by c = 1/y/JJpo and thus c ^ 1400 m.s~^. The friction term is 
given by the Manning-Strickler law (see 

Sf = K{S)U\U\ with K{S) ^ 



where S = S{X) is the surface area of the section r2(X) normal to the main pipe axis 
(see Fig. [T]for the notations). Kg is the coefficient of roughness and Rh{S) = S/Pm 
is the hydraulic radius where Pm is the perimeter of O. 

System ([9|) is integrated over the cross-section In the following, overlined letters 

represents the averaged quantities over Vt. For m G n = is the outward 

m 
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unit vector at the point m in the fi-plane and m stands for the vector torn (as 
displayed on FiG. d]). 

Following the work in [5], using the approximations pU ^ 'pU, pU"^ ^ 'pU and 
Lebesgue integral formulas, the mass conservation equations becomes: 

dtipS) + dxipq) = [ P (Udxm - V) .H ds, (10) 
Jdn ^ ' 

where q = S\J is the discharge of the flow and the velocity V = {V, WY in the 
(n, 5)-plane. 

The equation of the conservation of the momentum becomes 

dt{pq)+dx{^ + c^pS) = -gpS sine + c^p— 

- pSZ^igcosO) (11) 



+ / pU [Udxm - V j .n ds 

The integral terms appearing in (fTO]) and (fTTIl vanish, as the pipe is infinitely rigid, 
i.e. Q = Q{X) (see [5] for the dilatable case). It follows the non-penetration 
condition: 

U \ 

' .N = 0. 



Finally, omitting the overlined letters except Z, we obtain the equations for pres- 
surised flows under the form 

dt{pS)+ dxipq) = 

dtipq) +dx{^ + c^pS) = -pSgsinO - pSZ-^{gcos9) + c^p^ ^ 

where the quantity Z is the Z coordinate of the center of mass. 

Remark 2.3 In the case of a circular section pipe, we choose the plane curve 
{x,0,b{x)) as the mean axis and we get obviously Z = 0. 

pS 

Now, following [5], let us introduce the conservative variables A = — the equivalent 

Po 

wet area and the equivalent discharge Q = AU . Then dividing System (fT2l) by po 
we get: 

dt{A) + dx{Q) = 

dt{Q)+dx{^+c^A) = -gAsinO - Az4^{g cose)+ (1^) 

Remark 2.4 This choice of variables is motivated by the fact that this system is 
formally closed to the Shallow Water equations with topography source term in non 
uniform pipe. Indeed, the Shallow water equations for non uniform pipe reads [1]: 



( dtA + dxQ = 

+g cos 9I2 



dtQ + dx i^+gcoseh) = -gAsm9-A{h-h{A)/A)^{gcosi 
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where the terms glicosO, l2Cos9, {h — Ii{A)/A) are respectively the equivalent 
terms to c^A, c^A-^ln(5), Z in System ([13|). The quantities h, h, {h-Ii{A)/A) 

denotes respectively the classical term of hydrostatic pressure, the pressure source 
term induced by the change of geometry and the Z coordinate of the center of 
mass. Finally, the choice of these unknowns leads to a natural coupling between the 
pressurised and free surface model (called PFS-model presented by the authors in 



To close this section, let us give the classical properties of System (11311 : 

Theorem 2.1 (frictionless case) 

1. The system is stricly hyperbolic for A{t,X) > 0. 

2. For smooth solutions, the mean velocity U = Q/A satisfies 

dtU + dx(^ + HA/S) + + gz\=0 (14) 



where ^g{X) = / Z(^)— — cos6{^) for any arbitrary xq and Z the altitude 
Jxo dX 

1/2 

term defined by dxZ = sinO . The quantity — — h c^ln(^/5) + g^g + gZ is 
also called the total head. 
3. The still water steady states for U = is given by 

ln{A/S) + g^g + gZ = 0. (15) 
4- It admits a mathematical entropy 

E{A, Q) = ^ + c'^A ln{A/S) + gA^g + gAZ 
which satisfies the entropy inequality 

dtE + dxiiE + c^Ap) ^ 

Remark 2.5 

• If we consider the friction term, we have for smooth solutions: 

dtU + dx(^^+ HA/S) + g^e + gZ^ = -gK{S)U\U\ 

and the previous entropy equality reads 

dtE + dx{{E + c^A)U) = -gAK{S)U'^\U\ ^0 

• If we introduce Z the so-called pseudo altitude source term given by 

Z = Z + ^g ln(5) 

9 
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(where is defined in Theorem l2.1|) . we can rewrite System (fT3]l in the simpler 
form, close to the classical Shallow Water formulation: 



dtiA) + dxiQ) = 

dt{Q) + dx{^+p{X,A)] +gdxZ = 



(16) 



where p{X, A) = (? A. 
This reformulation allows us to perform an analysis close to the presented one by 
the autors in [9j in order to write the kinetic formulation. 



3 The kinetic model 

We present in this section the kinetic formulation (see e.g. [8]) for pressurised flows 
in water pipes modelized by System (fT6]l . To this end, we introduce a smooth real 
function x such that 

X(w^) = x(-w^) > 0, / x{w)dw = \, \ w^xiw)dw = 1 
Jk Jr 

and defines the Gibbs equilibrium as follows 



c \ c 



which represents the density of particles at time t, position x and the kinetic speed 
^. Then we get the following kinetic formulation: 

Theorem 3.1 The couple of functions {A,Q) is a strong solution of the Shallow 
Water-like system f7^) if and only if Ai satisfies the kinetic transport equation 

dtM + ^dxM-gdxZd^M = K{t,x,0 (17) 

for some collision kernel K{t,x,^) which admits vanishing moments up to order 1 
for a.e (t,x). 

Proof of Theorem 13.11 We get easily the above result since the following 
macro-microscopic relations holds 

A = [ MiOdC (18) 

JR 

Q = [ ^MiOdC (19) 

%- + ^A= f eM{i)di (20) 
^ Jr 

□ 

The reformulation of System (fT3]) and the above theorem has the advantage to give 
only one linear transport equation for M. which can be easily discretised (see for 
instance [HllTO]). Morever, the following results hold: 
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Theorem 3.2 Let us consider the minimization problem mm£{f) under the con- 
straints 

f>0, [ fiC)d^ = A, [ ^f{OdC = Q 
where the kinetic functional energy is defined by 

£{f) = I ^fiO + c^iOlogifiO) + c^f{0log{cV2^) + gZfiO 

A ( ^— U\ 

Then the minimum is attained by the function A4(t,x,^) = — x I I where 

c \ c J 

, . 1 (-w"' 



Morever, the minimal energy is 

£{M) = E{A,Q) = f- + c^AlnA + gAZ 
and M. satisfies the still water steady state equation for U = 0, that is, 

idxM- gdxZd^M = 0. 

Proof of Theorem 13.21 One may easily verify that / = 7W is a solution of 
the minimization problem. Under the hypothesis / > the functionnal £{f) is 
strictly convex which ensures the unicity of the minimum. Furthermore, by a direct 
computation, one has £{M) = E. 

The minimum M. of the functionnal £{f) satisfies the still water steady state 
for [/ = 0, 

idxM- gdxZd^M = 0. 

dxA ( ^\ A ( i\ 

Since dx-M. = x ( ~ ) ' ^5-^ = ~2X ( ~ ) ' denoting w = S,/c, we get 

~A 

wdxAxiw) - gdxZ—x'{w) = 0. 

c 

On the other hand, the still water steady state at macroscopic level is given by 

ln{A) + gZ = est, 

and so one has gdxZ = —c'^dxO-T^ A). Finally, we get the following ordinary differ- 
ential equation 

wxiw) +x{w) = 0. 

which gives the result. 

□ 



4 The kinetic scheme with pseudo-reflections 

This section is devoted to the construction of the numerical kinetic scheme and its 
properties. The numerical scheme is obtained by using a flux splitting method on the 
previous kinetic formulation (|17p . The source term due to the pseudo topography 
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dxZ is upwinded in a close manner described by Perthame et al. [9j using an 
energetic balance at the microscopic level. In the sequel, for the sake of simplicity, 
we consider the space domain infinite. 

Let us consider the discretization {niiji^z of the spatial domain with 



mi 



which are respectively the cell and mesh size for i G Z. Let At" = tn+i — tn, n & N 
be the timestep. 

Let = (^", Qf), = be respectively the approximation of the mean value 
of (A, Q) and the velocity U on rrii at time t^. 

Let A4"(^) = — ( ~ I be the approximation of the microscopic quantities 



and Zjlm-(X) be the piecewise constant representation of the pseudo-altitude Z. 
Then, integrating System (fT6ll over TTii X [t^, , we get: 

^r'=U--^{F-^,/,-Fl,,,) (21) 

where 



1+1/2 ^^n J V ^ ' ' 



^+l/2)dt] (22) 



are the interface fluxes with F{A, Q) = {Q, Q'^/A + (? A^- . 

1/2 



Now, it remains to define an approximation F^^,^ of the flux at the points XiJ^^ji- 



To this end, we use the kinetic formulation (fT7|) . 

Assume that the discrete macroscopic vector state is known at time tn- We 
consider the following problem 

/(t„,X,e) =A^(t„,X,0 (X,OGm, xM ^ ^ 

where A4(t„,X, ^) = A4"(^) in the cell mj. It is discretized as follows (since it is a 
linear transport equation) 



e z, Vn e N, /r+i(0 = Mm - i— {-^,"+1/2(0 - -Mti/2(0} (24) 



where denotes the interface density equilibrium (computed in section HH]) . 

Finally, we set 

Ur^ = J ( l^rHOdC (25) 



^ji+l / e _ Tjn+l 



and 



Remark 4.1 We can understand Equation (|23]) as follows: let us consider the 
following problem, 

f dtf + ^dxM - gdx{Z) S^TW = {t, X, ^ G [tn, Wi] x x R 

\ f{tn,X,0=M{tn,X,C) {X,C)emixR. ^ ' 
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Assuming that Ai{t,X,^) is known on x mj x M leads to the same dis- 

cretization ((241) of Equation (|23]) . Hence the numerical scheme ([241) avoids to com- 
pute explicitely the collision kernel K at the microscopic level. Indeed, substracting 
Equation (fT71) to Equation ([261) . we get: 

dt{M-fm = K{t,x,C). 
Then, integrating the previous identity in time t and ^ yields to: 

In other words, using the numerical scheme (l2^ and the macroscopic-microscopic 
relation ([25]) is a manner to perform all collisions at once and to recover exactly the 
macroscopic unknows {A,Q). 

Now to complete the numerical kinetic scheme, it remains to define the microscopic 
fluxes appearing in equation ([241) introduced by the choice of the constant 

piecewise representation of the pseudo-altitude term Z. 



4.1 Interface equilibrium densities 



To compute the interface equilibrium densities, we use the generalized characteristics 
method. Let s £ be a time variable and / the solution of the kinetic 

equation ([23l) . Let i £ Z, t £ and be respectively the kinetic speed 

of a particle at time t on each side of the interface Xj_|_i/2- The characteristic 
curves and X{s) of the kinetic transport equation ([^51) satisfies the following 
equations: 

^ = -gdJiXis)) 

(27) 

dX „ 
IF = 

where the final conditions are defined by 



X{t) 



X 



i+l/2 



(2^ 



for some constant ^ defined later. By a straightforward computation, we get the 
following mechanical conservation law: 



d_ 

ds 



+ gZis) =0 



(29) 



Since Z is a piecewise constant function, the solution H of the ordinary differential 
equation ([271) is a piecewise constant solution. So, we need to define an admissible 
jump condition to get only physical solutions of the problem ([27l) . Thanks to the 
relation ([29]) . we get the jump condition: 



2gZ 
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that is also 

5z 



2 2 



gAZi+y2 (30) 
where AZj_(_i/2 is such that 

Zj+i - Zi = AZj+i/2'^x,+i/2 

with (5a is the Dirac mass at point a. The quantity AZj_|_i/2 is the potential bareer. 
Next, solving System (f271) on mj x with the final conditions : 

= ^' , (31) 

X{t) = 

we get 

E{s) = 6 and X{s) = ^/(s - Wi) + ^*+i/2- (32) 



Reflection & Transmission 



1/2 



-< 



— Reflection 

Positive transmission 

~ Negative transmission 

□ Reflection zone 



Figure 2: The potential bareer: transmission and reflection of particle 
Top: the physical configuration 
Middle: the characteristic solution in {X, H)-plane 
Bottom: the characteristic solution in (X, t)-plane 

Due to the jump condition ([30]l and the sign of the kinetic speed, we distinguish 
three admissible cases as displayed on FiG. [2l 
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The case > corresponds to the positive transmission (this means that the 
particle comes from the left) and we deduce from Equalities (f32]l that the left 
microscopic flux -M^^-^^^^iO equal to M.^{^). 

The case ^/ < and — 2gAZi^i/2 < is the so-called reflection case. The 
condition — 2gAZij^i/2 < says simply that the slope of the X solution 



([321) cannot exceed y2(7AZj^x/2 (^s displayed on FiG. [2] (bottom)) and so 
the flux ■^^_^_i/2(0 given by Physically, since the particle with the 

kinetic speed under the previous kinetic condition, has not enough energy 
to overpass the bareer, it is reflected with the kinetic speed —^i. 

The last case is when ^ < and — 2gAZi_^_i/2 > 0. This case corresponds 
to the negative transmission: this means we take into account the particles 
coming from the right side with negative kinetic speed. Contrary to the re- 



flection case, the constraint on the X slope is limited by > —y 2gAZi_^_i/2 

and we get as solution A4"_)_x (^—\/ — ■ From a physical point 

of view, the observed particle at the left of the interface comes from the right 



side with a kinetic speed < where = ~y^f~ '^9^Zi+i/2j taking into 
account the gain or loss of potential energy through the bareer (as displayed 
on Fig. [l(bottom)). 

Finally, adding the previous results we obtain: 

positive transmission reflection 



negative transmission 



negative transmission reflection 



(33) 



positive transmission 

The microscopic flux at the right of the interface is obtained following a same 
approach. 



4.2 Numerical properties 

We present some numerical properties of the macroscopic scheme (f2T]) - (f22]) . namely 
the stability and the preservation of the still water steady state. The stability of 
the kinetic scheme depends on a kinetic CFL condition 

At" 

e < 1, ve 



maxj hi 



14 



and so, on the support of the maxwelhan function (e.g. we see that from the 
microscopic fluxes in Subsection 14. II) . The support of the maxwelhan function com- 
puted in Theorem 13.21 is not compact, then the stability condition cannot be sat- 
isfied. Therefore, in the sequel, we will consider the particular Gibbs equilibrium 

xiw) = ^] («^) introduced by the authors in p] and used in [6] in the 

case of pressurised flows in uniform closed pipe. 

Let us present the numerical properties of the scheme ((23]) - ([331) . 

Theorem 4.1 

1. Assuming the CFL condition 



max (ic/n +\/3c) < 1, 



maxjgz hi 

the numerical scheme (£^-(3^ keeps the wet equivalent area A positive. 
2. The still water steady state is preserved: 

f/" = 0, — ln(p") + Zi= est 
9 

Proof of Theorem 14.11 (It is similar to the one obtained in [S]) Let us suppose 
A" > for alH G Z and n S N. Let ^-t = max(0, be the positive or negative 

At" „ 

part of any real and a = — , Equation OSJ reads: 

maxj hi 

fr\o ^ (i-^ici)-Mr(o 



Since the support of the x function is compact, we get 

/f+i(e)>Oif <^/3c,VjEZ 

which implies |^| < +V^c. Using the CFL condition c7|^| < 1, we get the result. 
Morever, since /"^^ is a sum of positive term, we obtain /"^^ > 0, hence the wet 
equivalent area at time t"'~^^ is positive, i.e. 

To prove the second point, we distinguish cases ^ > and ^ < to show the equality 
■^1+1/2 = -^i-1/2- Using the jump condition (f30]) . we easily obtain /j""*"^ = Ai^ 
which gives the result. 

□ 
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Now let us also remark that the kinetic scheme (f2^ - ([33]l is wet equivalent area con- 
servative . Indeed, let us denote the first component of the discrete fluxes (^yl)^i/2' 

iFA)ti/2--= [ iM%,,,{Odi 

An easy computation, using the change of variables nP' = S^"^ — 2glS.Zi^xl2 in the 
interface densities formulas defining the kinetic fluxes ^Ho'^s us to show 

that: 

5 Numerical Validation 

The validation is performed in the case of a soft and sharp water hammer in an 
uniform pipe. Then we compare the results to the ones provided by an industrial 
code used at EDF-CIH (Prance) (see [E]), which solves the Allievi equation by the 
method of characteristics. The validation in non uniform pipes is performed in the 
case of an immediate flow shut down in a quasi-frictionless cone-shaped pipe. The 
results are then compared to the equivalent pipe method [T]. 

5.1 The uniform case 

We present now numerical results of a water hammer test. The pipe of circular cross- 
section of 2 m^ and thickness 20 cm is 2000 m long. The altitude of the upstream 
end of the pipe is 250 m and the slope is 5°. The Young modulus is 23 10^ Pa since 
the pipe is supposed to be built in concrete. The total upstream head is 300 m. 
The initial downstream discharge is 10 m^/s and we cut the flow in 10 seconds for 
the flrst test case and in 5 seconds for the other. 

We present a validation of the proposed scheme by comparing numerical results of 
the proposed model solved by the kinetic scheme with the ones obtained by solving 
Allievi equations by the method of characteristics with the so-called belier code: 
an industrial code used by the engineers of the Center in Hydraulics Engineering of 
Electricite De Prance (EDP) [E]. 

A simulation of the water hammer test was done for a CPL coefficient equal to 0.8 
and a spatial discretisation of 1000 mesh points. In the figures PiG. [3] and PiG. 
m we present a comparison between the results obtained by our kinetic scheme 
and the ones obtained by the "belier" code: the behaviour of the discharge at the 
middle of the pipe. One can observe that the results for the proposed model are 
in very good agreement with the solution of Allievi equations. A little smoothing 
effect and absorption may be probably due to the first order discretisation type. 
A second order scheme may be implemented naturally and will produce a better 
approximation. 
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Figure 3: Comparison between the kinetic scheme and the industrial code belier 
First case: discharge at the middle of the pipe 
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Figure 4: Comparison between the kinetic scheme and the industrial code belier 
Second case: discharge at the middle of the pipe 
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5.2 The case of non uniform circular pipe 



We present a test of the proposed kinetic scheme in the case of a contracting or 
expanding circular pipes of length L = 1000 m. The downstream radius is kept 
constant, equal to i22 = Iw- and the upstream radius varies from Ri = Im to Am 
by steps of 0.25 m. The others paramaters are = 300 mesh points, Ks = 9000 
(this means that the wall of the pipe is very smooth), CFL= 0.8. The upstream 
discharge before the shut-down (1.5 seconds) is fixed to lOm^.s"^ while the up- 
stream condition is a constant total head. We assume also that the pipe is rigid. 
Then for each value of the radius -Ri, we compute the water hammer pressure rise 
at the position x = 96 m of the pipe and we compare it to the one obtained by the 
equivalent pipe method (see [T]). The results are presented in FiG. [S]and show a 
very good agreement. 

We point out that the behaviour of the solutions corresponding to the equivalent 

pipe method and our method are different: this is due to the dynamic treatment 
d 111 S 

of the term (? — -— — related to the variable section which is not present in the 

dX _ _ _ 

equivalent pipe method: see FiG. El FiG. d FiG. [8]: 
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Figure 5: Comparison in the prediction of pressure rises in cone-shaped pipes between 
the present method and the equivalent pipe method 
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Discharge at x = 96 m; Rl = 1.25, R2 = 1 



Piezometric line at x = 96 m; Rl = 1.25, R2 = 1 
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Figure 6: Discharge (left) and piezometric line (right) for Ri = 1.25 m 
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Figure 7: Discharge (left) and piezometric line (right) for Ri = 1.5 m 
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Figure 8: Discharge (left) and piezometric line (right) for Ri = Am 
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